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1 Introduction 



Multi-jet final states play an important role in the study of high-energy collisions. They provide 
in fact interesting signatures for several phenomena, both within the Standard Model (e.g. top- 
pair production), and beyond it (e.g. multi-jet decays of supersymmetric particles such as gluinos 
and squarks). The accurate determination of the properties of these phenomena requires a good 
understanding of the properties of the usually large multi-jet QCD backgrounds, which can distort 
the shapes of signal distributions and affect the measurement of quantities such as resonances' 
masses. In the past few years, several theoretical developments have allowed the calculation of very 
complicated multi-parton processes in QCD (for a review, see 0) 0. For example, the complete 
set of leading-order (LO) background processes to the production and decay of top-quark pairs in 
hadronic collisions is known for both the fully hadronic decays 0, and the ez/+4-jet decays [[|. 

While significant progress has been made in the field of one-loop corrections (for a review 
see ||), quantitative studies of multi-jet processes (with n > 3) can only be done today using 
tree-level results, which are the subject of the present work. There are several reasons for wanting 
to improve the tools currently available to perform these calculations. 

1. First of all, interesting final states with larger jet multiplicities will become available with 
the next generation of colliders (LHC and NLC). This will hold both for standard QCD 
processes, where the huge available phase-space will allow production of many high-.&r 
jets, and for potential signals of new physics (a good example being cascade decays of heavy 
squarks or gluinos with i?-parity non-conservation, where final states with over 10 jets would 
be a typical signature). This calls for improved algorithms, to allow the calculation of cross- 
sections for complicated processes within reasonable amounts of computer time. 

2. Secondly, one would like to be able to complement the calculation of parton- level matrix 
elements with the evaluation of the full hadronic structure of the final state. This is a 
key ingredient for a satisfactory study of both signal and background components and for 
a complete comparison between theory and data. It involves the consistent merging of 
the matrix-element computation with the parton-shower evolution, a problem which has not 
been considered in the development of previous tools for the evaluation of multi-jet processes. 

While relations are known || which allow to systematically evaluate high-order tree-level 
processes in a recursive fashion, the complexity of the algorithm grows very quickly and makes 
progress beyond the processes listed above very hard. Several approximations have therefore been 
introduced || [7], to evaluate with rather low computing time and acceptable accuracy cross- 
sections for many partons in the final state. Nevertheless, the knowledge of the exact parton level 
matrix elements will always be needed, in order to assess the reliability of the approximations 
used. 

In this work we present an approach to the calculation of tree-level matrix elements for multi- 
parton final states which addresses both of the above problems, and therefore improves on the 
currently available algorithms. The key element of this proposal is the use of the algorithm ALPHA, 
developed by two of us H for the evaluation of arbitrary multi-parton matrix elements. This 
algorithm determines the matrix elements from a (numerical) Legendre transform of the effective 

3 We shall include under the label QCD also all processes with the emission of electro-weak gauge bosons, such 
as associated production of W's and jets. 
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Process 


n = 7 


n = 8 


n = 9 


n = 10 


9 9 ~>ng 


559,405 


10,525,900 


224,449,225 


5,348,843,500 


qq->ng 


231,280 


4,016,775 


79,603,720 


1,773,172,275 



Table 1: Number of Feynman diagrams corresponding to amplitudes with different num- 
bers of quarks and gluons. 



action, using a recursive procedure which does not make explicit use of Feynman diagrams. The 
algorithm has a complexity growing like a power in the number of particles, compared to the 
factorial-like growth that one expects from naive diagram counting. This is a necessary feature of 
any attempt to evaluate matrix elements for processes with large numbers of external particles, 
since the number of Feynman diagrams grows very quickly beyond any reasonable value. For 
example, the number of tree-level Feynamn diagrams corresponding to a process with n g gluons 
and n q pairs of quarks (with n q = 0, 1), is given by the following formula 



diag 



y7r + x y 3 7r + z2 yjr) x ^ y ^ z ^ 

. ox ay oz I 



(1) 

-y=z=l 



A similar formula can be constructed for processes with more than 1 quark pair. The number 
of diagrams relevant for some of the examples considered in this paper are collected in Table 1. 
These numbers clearly illustrate the problems encountered when trying to evaluate the amplitudes 
by calculating each individual diagram, whether by algebraic or by numerical means. 

ALPHA has been shown to operate very successfully in the case of purely electroweak pro- 
cesses ||, and will be reviewed here in Section 2. Its application to the case of QCD, although 
straightforward, requires some care due to the rapid increase of the number of colour configurations 
with the number of external coloured particles. Furthermore, the choice of how to organise the 
sum over colour structures has important implications for the possibility to merge the evaluation 
of a given matrix element with the successive parton-shower QCD evolution. Several options are 
available, in principle, to deal with the summation over all possible colour configurations. They 
will be discussed in Section 3, where our strategy for an efficient generation of Monte-Carlo events, 
combined with the possibility to generate colour configurations suitable for a parton-shower evo- 
lution, will be presented. Some numerical examples of applications of this tecnique to the case 
of 2 — > 8 parton processes in hadronic collisions are given in Section 4. A more complete study, 
including the treatment of associated production of jets and electroweak gauge bosons and the de- 
scription of a complete event generator interfaced to a parton-shower program such as HERWIG [ 110] , 
will be presented in the future. 



The formula can be easily derived as follows. For a given diagram consider the triplet (p, q,g), where p is the 
number of 3-gluon vertices, q is the total number of external quark legs and internal quark propagators, and g is 
the total number of gluonic external legs and internal propagators. The number of diagrams obtained by attaching 
an additional gluon in all allowed ways is given by: p diagrams of type (p — l,q,g + 1), plus q diagrams of type 
(p, q + 1,5+ 1) plus g diagrams of type (p + l,q,g + 2). These triplets can be obtained by applying the operator 
yd I dx + xy 3 d / dy + z 2 yd / dz to the function x p y 9 z q . Their multiplicities are given by the coefficients of the relative 
monomials, and the total number of generated diagrams is extracted by setting x = y = z= lin the polynomial 
thus obtained. Repeated iteration of this operator on the three-point diagrams ggg and qqg, represented by the 
monomials xy 3 and y z 2 , gives the desired result. 
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Independently of our efforts, work by Draggiotis, Kleiss and Papadopoulos |TT] has recently 
addressed the problem of the efficient generation of multi-parton QCD final states using the 
ALPHA algorithm. In this work they turn the summation over colours into an integration over a 
continuous set of variables, slightly gaining in computational accuracy over the more standard 
colour-summation technique presented in our work. The integration technique, however, does not 
lend itself in an obvious way to the efficient merging of the parton-level calculation with the parton- 
shower evolution. In either case, it is quite clear that improvements in the numerical efficiency 
of these calculations will be possible, and work should be devoted in the future to find the best 
compromise among all different requirements to be met for a faithful and efficient representation 
of these highly complex processes. 



2 The ALPHA algorithm 

In reference ||, a new approach to the computation of tree level scattering amplitudes was in- 
troduced. This approach, based on the numerical Legendre transform of the effective action, is 
particularly useful for the automatic calculation of multi-particle processes. This technique was 
implemented in a Fortran code || which has been succesfully used to study several intricated 
electroweak processes ||. For the sake of completeness, we review in this Section the ALPHA al- 
gorithm; the interested reader can find a more detailed discussion in the original paper || , which 
includes an explicit analytic example for the A0 3 theory. 

Let T be the one-particle-irreducible generator of the Green functions for a given theory. Then 
the computation of the S-matrix requires the evaluation of the Legendre transform, Z, of V: 

z{F) = -r(<n + j a ( X )r(x) (2) 

where <p a are the classical fields defined as the solutions of 

ja = ?jr> (3) 

and the J a play the role of classical sources. In general the Lagrangian contains several fields, with 
different spin and internal quantum numbers. It may also contain interactions with an arbitrary 
number of fields; for our purposes it is better to rearrange the Lagrangian into an equivalent form 
which includes only trilinear interactions: this can be achieved by introducing a proper set of 
auxiliary fields^]. The Lagrangian in momentum space is 



£ = \jd 4 p d 4 q 6(p + q) fl af3 (p 2 ) <j> a (p) <fP(q) + 



+ - J d 4 p d 4 q d 4 k 5(p + q + k) <p a (p) <f>P(q) ^(k) O a ^(p, q, k) , 

where the greek indices are a compact notation for Lorentz indices, internal symmetry indices, 
flavour etc. and where d 4 p = d 4 p/{2ir) 4 and 8 = (27r) 4 <5 4 . U al3 (p 2 ) is the inverse propagator 
and O a P~> [p,q,k) is a generic function of the momenta. In the case of translationally invariant 

5 For example a term of the type A</> 4 can be replaced by the equivalent form X(f> 2 X — 1/4X 2 , where X is an 
auxiliary field. The equation of motion for the auxiliary field, X — 2A^> 2 , makes this equivalence manifest. 



3 



local interactions O is a polynomial in the momenta. To obtain the connected Green function 
n ) from the above Lagrangian we introduce the classical sources 

n 

J a {q)=T.^-Pi) ( 4 ) 

i=l 

where the a" carry the same quantum numbers as the source J a . With this choice of the J a the 
amplitude A is given byf| 

dZ , 

A ~ da?...daZ la "=°- aZ=0 (5) 



The equations of motion of eq. (Q) can be written as 



J x (q) - - / d 4 p d 4 k 5{q + k + p) <f>P(k) <f{p) X ^(q, k,p) 



(6) 



It is clear that this equation for <p(q) can be solved perturbatively with respect to the interaction 
Qap-y ^ or? equivalently, with respect to the af) and the (t + l)-th order of this perturbative series 
is obtained by inserting the expansion of the 4> a (q) up to the t-th order in the right-hand side of 
eq. (H). To recover the functional derivative of eq. (|5|) and avoid unnecessary computations it is 
useful to introduce the prescription to drop out terms which contain powers of the a" larger than 
one, at any iteration step. With this prescription, using the initial condition (|]) and the recursive 
relation (|B]), we can prove by induction that the solution (f) a (q) is of the form 

2 n — 2 

with 

P 3 = c)pi 4 = 0,1. (8) 

Here an important point should be noticed: the Lagrangian £ of eq. (f|) is reduced to a simple 
function of a finite number of 6" variables, by means of eq. (d) with the constraint given in eq. (§). 
Namely, plugging the explicit expressions of eq. (|7]) into the Lagrangian of eq. (Q) we obtain 

£ = l - J d 4 p d 4 q S(p + q) n Q/3 (P/) 5(p - Pj) 5(q - P r ) b« b% + 

+ \\d 4 p d 4 q d 4 k 5(p + q + k) O a ^(p, q, k) 5(p - Pj) 5(q - P r ) 5(k - P t ) b« ^ bj (9) 

We define the matrices 

Af = J d 4 p d 4 q S(p + q) fH*{Ff) 5(p - Pj) S(q - P) (10) 

and 

D^Z = J d 4 p d 4 q d 4 k 8{p + q+ k) O a ^(p, q, k) 5{p - P 3 ) 5{q - P) 8{k - P m ). (11) 
6 For the sake of simplicity, we omit from this expression the explicit truncation of the external propagators. 
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After integration they become: 



O a ^{P 3) P u P m ) if P j + ft + P m = 



if P j +P l + P m ^0 

fir^(if) if p d + Pi = o 

A? = ' (13) 

{ if Pj + P^O. 

We recall that, in the above expressions, O alBl is the interaction in momentum space, U a/3 is the 
inverse propagator and the quantities Dji m and Aji satisfy the four-momentum conservation by 
construction. 

The function Z then becomes 

Z = <W - i b« bf Af - i 6? 6f % D/J. (14) 

Equation (|6]) provides us with a set of iterative relations for the 6": after the substitution in 
eq. (|7|), and performing the integrals over the momenta, we are left with a relation which gives us 
each bj at the order t (6j,t) in the interaction coefficient O, in terms of the bj >r (with r < t)[] 

Namely, we have |§ 

6? = aj j = l,n 

= j > n n— number of external particles 

bh = -^(A-^K^E^to (15) 



bj,t — 9 (A 1 )jlL^m,k,l X/ ^fc,r ^, 
Z r+s= t_l 



where the condition / 7^ derives from the constraint in eq. (§). The full scattering amplitude is 
recovered by plugging 6" = ^ &j,t m ec l- (0) an d keeping only terms which are proportional to 
a^af ■ ■ ■ o-Z (the only ones that contribute to the limit a — > 0) 

1 1 

/I - _Z V /) a A Q/3 ^ — - V D^ 6 y 3 K 1 h s + h a n 13 T\ alS flfTl 

Z s+r=n-2 s+r+t=n-3 

Notice that in eqs. (|i~4"|) -( fL6|) repeated indices are summed over, in particular repeated latin indices 
are summed from 1 to 2 n — 2, namely over the non-zero momenta contributing to eq. (0). 

Equation (|T^) can be further simplified using the equations of motion. The solution bj to 
dZ/dbj = 0, has been found for any value of the expansion parameter, thus the minimization is 
satisfied order by order in the perturbative expansion, i.e. dZ/dbjj = 0. For simplicity, let us 
discuss the case of an odd number of particles n. It is easy to check by inspection that each term 
in eq. (|i~6l) can contain at most one bj t t with t > n/2 — 1. Therefore we can reorganize eq. (|i~6f) by 

7 In the following j, k, I, m will always label quantities that are in one to one correspondence with the momenta 
Pj in eq. (^), whereas r, s, t will denote the order of the perturbative expansion. 
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collecting each hj t with t > n/2 — 1. The coefficients of these terms cannot depend on other bjj' 
with t' > n/2 — 1 (they are only polynomial in bjj with £ < n/2 — 1) and therefore we can simply 
drop such terms from eq. (|T6|). The above simplification also applies to the case of even n, and it 
amounts to keeping only the trilinear terms in eq. flTBP and limiting the summation to a subset of 
the bj ||. Then, an additional restriction on the quantities c*- in eq. (|8|) is obtained: 

Ai,-,Pn = -I E E 6f,r &Z, , (i7) 

with: 



6 

i,fc,ZeP«+r-+t=n-3 



E*4<n/2 

j EV it { or (18) 

Ei 4 = n/2, with c) = 1 . 

A final remark is in order here. In eqs. (jl5[) one has to take properly into account Bose/Fermi 
statistics. Formally this can be achieved by introducing a set of Grassman variables 6j, so that 
£j£k + £k£j = 0, and setting u(p) — > eu(p) for the fermion sources, where u(p) is a vector of ordinary 
numbers. In practice this means that each term in the sums of eqs. ( |T5| ) enters with a relative 
sign, depending on the order of the bj. 

The advantage of the iterative eqs. flT5| ) is that they can now be easily implemented in a 
Fortran code. Note that the recursive relations have been obtained for a completely generic 
Lagrangian, without using any specific property of the interaction (e.g. identities on the structure 
constants of SU(3), or the Lorentz structure of the interaction etc.). With this algorithm we can 
thus compute the scattering amplitude for any physical initial and final states. In particular, 
any colour structure can be assigned to the external legs (while for example the Berends-Giele's 
recursive relations Jl5]| have only been derived for colour-ordered amplitudes) and weak bosons (as 



well as other particles beyond the Standard Model) can be incorporated. Finally, the algorithm 
has an exponential growth of the CPU time with the number of external particles (instead of a 
factorial growth) and it can take into account the particle masses without increasing the computing 
time. 

As a final remark, we point out that the algorithm presented above can be further optimised 
in the specific case of the QCD Lagrangian. In fact the pure Yang-Mills Lagrangian 

Cym = -^F^ V F* U (19) 

can be rewritten as 

C-ym = —^B a iW B'/ w + -B° U F° V (20) 

This form of the Lagrangian has the virtue that a single interaction term of the form BAA is left, 
instead of the two interaction terms which are present in eq. (|T9|) . This saves CPU time when 
performing the iteration step given in eq. (15). A further reduction in the algorithm complexity 
can be obtained by using the Coulomb gauge Aq = 0, which reduces the number of components 
of the field A ^ from four to three and those of the field from six to three. 



3 Summing over colours 

The proper bookkeeping of the colour structure of QCD processes has been one of the main 
ingredients in the simplification of the calculations of multi-parton processes in QCD which took 
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place in the past years. It has been shown |12|| that by expanding the matrix element for a given 
QCD process in a particular colour basis, the coefficients of each colour structure in this basis enjoy 
important properties which make their calculation much simpler. As an example, the scattering 
amplitude for n gluons with momenta pf, helicities ef and colours a,i (with i = 1, . . . , n), can be 
written as: 

M({ Pi }, {e 4 }, {a,}) = £ tr(A ai A" 2 . . . \ a ") A{1, 2,...,n). (21) 

P(2,3,...,n) 

The sum extends over all permutations P of (2, 3, . . . , n), and the functions A(l, 2, . . . , n) (known 
as dual or colour- ordered amplitudes) are gauge-invariant, cyclically-symmetric functions of the 
gluons' momenta and helicities^]. All of the colour-dependence is absorbed in the trace coefficients, 
and the dual amplitudes are colour-independent. In these special bases some particular helicity 
amplitudes |14[ have very simple analitical expressions [|1^, [5], regardless of the number of external 
partons, and all amplitudes obey recursive relations || which make it possible to numerically 
evaluate them systematically for arbitrarily complex processes |0| [15], Thanks to additional 
properties of these functions, only (n — 2)! of them are truly independent, and the remaining ones 
can be obtained from specific linear combinations |]13|| . 

Furthermore, when summing over colours the amplitude squared, different orderings of dual 
amplitudes are orthogonal at the leading order in 1/N 2 : 



E MM, {*}, M)| 2 = n-\n 2 - i) E 

col's P(2,3,...,n) 



|A(l,2,...,n)| 2 + ^(interf.; 



(22) 



It is therefore possible to achieve an accuracy to leading-order in 1/N 2 by neglecting the evaluation 
of the subleading interferences, reducing significantly the complexity of the numerical evaluations. 

Similar expansions hold for processes involving one []T7| , H or more Jl8| quark pairs. In the 
case of amplitudes with one quark pair, for example, one has the following expansion: 

M{q a ,qp,{pi},{ei},{<k}) = E (A fll A fl2 . . . A a ")^A(g, g, 1, 2, . . . , n) (23) 

P(l,2,...,n) 

where the gauge-invariant functions A(q, q, 1, 2, . . . , n) obey remarkable properties similar to those 
of the gluonic amplitudes. 

Dual amplitudes can be easily evaluated using the ALPHA algorithm. Since the dual amplitudes 
A are independent of the number of colours, they can be calculated exactly by taking iV sufficiently 
large. Considering for example the case of an n-gluon amplitude, we can choose N > n and select 
the following set of A matrices to represent the gluon colours ai, . . . ,a n : 

(A ai ) jfc = W (« = 1, . . ■ , n - 1) , (X an ) jk = ^S n>j 5 k>1 (24) 

With this colour choice the dual amplitude corresponding to the permutation (1,2, ... ,n) is pro- 
portional to the full amplitude, as the only non-vanishing colour factor in eq. ([H]) is 

tr(A ai A a2 ...A ffl ") = J^. (25) 



8 For simplicity we will just use the indices i = 1, . . . , n, as opposed to using the full symbols pi and e», to specify 
the relevant permutation of momenta and helicities. 
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Figure 1: Colour structure of the n-gluon amplitude in the large- N limit. 



We explicitly calculated multi-gluon dual amplitudes using the large-iV colour-matrices given 
in eq. fl24|). We verified the correctness of the calculation for n up to 11 by comparing the 



results for maximally helicity violating (MHV) amplitudes |TJ| (e.g. g + g + — > g + ■ ■ ■ g + ) with the 



analytic expressions known exactly for arbitrary n [12|, ||. The input of the numerical evaluation 
of the matrix element is a string containing the total number of gluons, their helicity state, and 
their momenta. From these data, the amplitude is evaluated automatically. Since the evaluation 
performed with the ALPHA algorithm does not treat the case of MHV amplitudes differently than 
any other helicity combination, and since we checked the numerical calculation using several 
different gauges (which allowed us to ensure that no accidental cancellation of individual diagrams 
takes place), we are confident that our code correctly evaluates the dual amplitudes for arbitrary 
helicity configurations. We extended this test to processes with one and two qq pairs, where 
analytic expressions are known for similar MHV amplitudes and found perfect agreement. 

The use of dual amplitudes for the evaluation of multi-parton processes has one valuable 
feature, and one serious drawback. The valuable feature is the fact that dual amplitudes admit 
a simple physical interpretation, which as we shall show makes them the required starting point 
for the parton-shower evolution of the parton-level process. The serious drawback is that their 
number grows factorially with n, so that while each individual dual amplitude can be easily and 
quickly calculated, using the technique given above, the number of dual-amplitude calculations 
which are necessary to get the full matrix element squared and summed over colours becomes soon 
too large to be practical. In the rest of this Section we will first explain in more detail the role of 
dual amplitudes for the parton-shower evolution of the hard process, and then present a way to 
bypass the rapid growth of their multiplicity. 

Dual amplitudes correspond to planar amplitudes in the N — > oo limit of QCD. At large N, 
the colour structure given by the assignement in eq. (H) corresponds to having the colour flowing 
from gluon 1 to gluon 2, from gluon 2 to gluon 3, and so on, until the colour of the last gluon 
flows back to gluon 1 (see fig. 1). The identification of a specific colour flow makes it possible to 
incorporate soft-gluon emission corrections to the hard process. In fact the soft-gluon emission 
probability from a planar amplitude is given, in the large- N limit, by the incoherent sum over the 
emission probabilities from each individual colour-string 
amplitude, for example, we get: 



IS, 12, [L5[. In the case of a multi-gluon 



Y d \M{p x ,...,VnM 



2 fc^O 2 



9 2 n E E 



{Pi ■ V 



i+i, 



col's 



\A(l,2,...,n)\ 2 + 0(l/N 2 ) 



(26) 

where n + 1 and 1 are identified in the above equation. Inclusion of soft-gluon virtual corrections 
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factorises in a similar fashion, and Sudakov form factors for the soft-gluon emission probability 
can be defined, to next-to-leading logarithmic accuracy, to describe the parton-shower evolution 
from any such colour-ordered process [^, [TO]. Perfectly similar equations can be written in the 
case of processes involving quark pairs. 

The prescription to correctly generate the parton-shower associated to a given event in the 
large- N limit is therefore the following: 



1. Calculate the (n — 1)! dual amplitudes corresponding to all possible planar colour configu- 
rations. 

2. Extract the most likely colour configuration for this event on a statistical basis, according to 
the relative contribution of the single configurations to the total event weight (]. Since each 
dual amplitude is gauge invariant, the choice of colour-configurations is also a gauge- invariant 
operation. 

3. Develop the parton shower out of each initial and final-state parton, starting from the 
selected colour configuration. This step can be carried out by feeding the generated event 
to a Monte-Carlo program such as HERWIG, which is precisely designed to turn partons into 
jets starting from an assigned colour-ordered configuration. 



Notice that, if the dual amplitudes are evaluated for a specific helicity configuration, HERWIG will 



also include spin-correlation effects in the evolution of the parton shower [^, [10]. 

With the physical value of iV = 3, dual amplitudes corresponding to different permutations 
will however interfere with each other, as shown in eq. (|22|). This interference is suppressed by 
powers of 1/N 2 Jr7[] , as well as by dynamical factors: for example, the behaviour of interference 
terms is less singular near collinear or soft momentum configurations than the leading terms in 
N. Within the 1/N 2 approximation which is employed in the description of coherence effects in 
the shower evolution |TIJ, it is therefore consistent to neglect the interferences between different 



dual amplitudes for the selection of the colour structure to be assigned to a given event. After 
calculating the total weight of a given event, accurate to all orders in 1/N, we can thus still use 
the procedure described in point 2 above to assign a definite colour configuration to the event 
itself. 

As a result, use of the dual-amplitude representation of a multi-gluon amplitude allows to 
accurately describe not only the large-angle correlations in multi-jet final states, but also the full 
shower evolution of the initial- and final-state partons with the same accuracy available in HERWIG 
for the description of 2-jet final states. 

In the N — > oo limit the choice of the dual amplitude basis has also one important advantage. 
Since interferences between different colour structures vanish in this limit, one can perform the sum 
over colours by Monte-Carlo methods. Rather than evaluating the full matrix element squared, 
summed over all colour structures, one can randomly select a dual colour structure on an event- 
by-event basis, and just evaluate the corresponding contribution to the amplitude squared. An 
overall multiplicative coefficient proportional to the number of dual colour configurations provides 
the correct normalization. Assuming that, on average, all colour configurations contribute the 
same amount to the cross-section, this approach is numerically more efficient than summing each 



9 Dcfining Wi = \Ai\ 2 for each colour flow i, and Wi = J2k=i » Wk l J2k=i n w k> tne J~th colour structure will 
be selected if Wj-i < £ < Wj, for a random number £ uniformly distributed over the interval [0, 1]. 
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event over all colours. Furthermore, one could optimise the selection of colour configurations, to 
adapt it to possible differences in their individual overall contributions. This is similar to what is 
usually done to perform the sum over quark and gluon helicities. 

At finite N this procedure is not applicable anymore, as interferences between various colour 
structures do not vanish. At the same time, the matrix describing all possible colour interferences 
has a size growing like [(n — l)!] 2 , which makes its storing and access highly inefficient. 

We propose to solve this problem as follows. First of all we choose the following orthonormal 
basis for the Gell-Mann A matrices: 



A 1 



1 

Til 

' - 

/ 

1 

V o o o 





A 2 = -4= 



1 

, / 1 



-1 

V o o 









A 3 = — = 



A 7 




In this basis, only a fraction of all possible 8 n colour assignements gives rise to a non-zero am- 
plitude. For each event, we randomly select a non- vanishing colour assignement for the exernal 
gluons, and evaluate the amplitude M. The weight of the event is proportional to |M| 2 , multiplied 
by the number of non-zero colour configurations. This is all we need if we are not interested in 
evolving the event with the parton shower. If instead we want to generate the parton shower, we 
first decide, with standard unweighting techniques, whether to accept the event. If the event is 
accepted, we list all dual amplitudes contributing to the chosen colour configuration according to 
eq. ([H]) and, among these dual amplitudes, we randomly select a colour flow on the basis of their 
relative weight^. 

In a 6-gluon amplitude, for example, a possible non-zero colour assignement is given by 
(2, 7, 5, 6, 1, 3). Up to cyclic permutations, only three orderings of the colour indices give rise to 
non- vanishing traces: tr(A 2 A 7 A 5 A 6 A 1 A 3 ), tr(A 2 A 6 A 1 A 5 A 7 A 3 ) and tr(A 2 A 7 A 3 A 1 A 5 A 6 ). Therefore 
only three dual amplitudes contribute to the full amplitude: A(2, 7, 5, 6, 1, 3), A(2, 6, 1, 5, 7, 3) and 
A(2, 7, 3, 1, 5, 6). The colour ordering to be specified for the coherent parton-shower evolution can 
be selected by comparing the size of the squares of tr(A 2 A 12 . . . A* 6 ) A{2, . . . , z 6 ) for the three 
contributing permutations (z 2 , . . . , ie) of the colour indices. 

Since the average number of dual amplitudes involved in the evaluation of a single element of 
the orthogonal basis is smaller than (n — 1)!, the complexity of the procedure grows more slowly 
than for the calculation done using directly the dual basis. The distribution of the number of dual 
amplitudes contributing to all possible colour assignements in the orthogonal basis for n = 8, 9 
and 10 gluons is shown in fig. 2. Furthermore, considering that only unweighted events are usually 
evolved by the parton shower, and that unweighted events are a small fraction of all generated 
parton-level events, the decomposition in terms of dual amplitudes only needs to be performed 
for a small fraction of the generated configurations. 



10 We explicitly checked the numerical implementation of this algorithm by comparing our results for the colour- 
summed squared amplitudes of the gg — > ng and qq — > ng processes (for n = 4,5) with the known results obtained 
in ref. |l3), as implemented in the NJETS code |||. The agreement is at the level of machine precision. 
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Figure 2: Distribution of the number of dual amplitudes contributing to all possible 
colour assignements in the orthogonal basis for n = 8, 9 and 10 gluons. 

The results we present in the following are all relative to parton-level results, and therefore 
only the Monte-Carlo summation over orthogonal colour configurations is considered. 



4 Results 



As an example of our technique, we present here results for the following two parton-level processes: 



9 9 
q q 



8 9 
8g. 



(27) 



For comparisons, we also computed the above reactions in the simple approximation first sug- 
gested by Kunszt and Stirling 0. This approximation (hereafter referred to as SPHEL) consists 
in assuming that the average value of MHV amplitudes is equal to the average value of all other 
non-zero amplitudes. In the case of n-gluon amplitudes this amounts to estimating the sum over 
all helicity configurations using the relation: 

E \M(gg^(n-2)g)\ 2 = I + > £ \M(gg^(n-2)g)\ 2 (28) 

hel's ) MHV 

where the sum on the right-hand side runs over all MHV amplitues (e.g. +H ►+...+, H > 

— }-...+, etc.). Their value is known exactly for all n at the leading-order in 1/N [p^j : 

n-2 



X; \M{gg^{n-2)g) 



9 2 N' 



MHV 



h3 
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Process 


hxact (pb) 


SPHEL (pb) 


9 9 ^8 g 


0.719 ± 0.019 


1.53 ± 0.03 


q q -> 8 g 


(0.901 ± 0.009)xl0- 3 


(1.042 ± 0.008 )xl0~ 3 



Table 2: Partonic cross-sections in pb at \/§ = 1500 GeV, with the cuts described in the 
text. 

1 



y 

P( W) (Pi ' P 2 ) (Pa • Ps) • • • (Pn • Pi) 



(29) 



where the sum runs over all permutations P(l, . . . , n — 1) of the indices (1, . . . , n — 1). Similarly, 
in the case of ggg . . . g amplitudes, the SPHEL approximation amounts to assuming: 



E \M 

hel's 



in— 1 



n ^)r 



n 



E \M(qq^ng)f 



MHV 



where the sum on the right-hand side runs over all MHV amplitues (e.g. +- 
value is known exactly for all n at the leading-order in 1/N ||17|| : 



(30) 

+ ...+), whose 



( a 2 N\ " 

E \M(qq^ng)\ = 4 (iV 2 - 1) £ [(p, • p,) 3 (^- • Pi ) + (p q ■ Pi )(p g ■ Pi f 

MHV V A ) i=l,n 

1 1 

Pi ' P? p(w) ^ 9 ' ' P 2 ) ' ■ ■ (p» ' p<?) ' 

where the sum runs over all permutations P(l, . . . , n) of the indices (1, . . . , n). 

The kinematic configuration and the cut values used in our numerical examples are as follows: 

VI = 1500 GeV , p Ti > 60 GeV 



\r]i\<2, Ai^>0.7. 



(32) 



These values, and the choice of a fixed strong coupling a s = 0.12, only serve for illustrative pur- 
poses. A more complete phenomenological analysis of production cross-sections and a comparison 
of exact and approximate expressions will be presented elsewhere. 

The integration over the phase space allowed by the cuts was performed by Monte Carlo using 
both RAMB0 |21| , a flat phase space generator, and a multichannel approach |22| . Both integration 
strategies gave compatible and comparable efficiencies. The sum over helicity configurations was 
performed via a flat Monte Carlo generation. No attempt has been made to optimise this step. 
Some of the details of the numerical performance of the algorithm are given at the end of the 
Section. 

In Table 2 we present our Monte Carlo estimate for the partonic cross-sections, together with 
the values obtained by using SPHEL. SPHEL overestimates by about a factor 2 the exact cross-section 
for the process g g — > 8 g , while it is accurate at the 10% level for q q — ► 8 g. 

Notice the large ratio of the gg-initiated amplitude, relative to the qq one. In the case of 2 — > 2 
processes, this ratio is of order 30 for a reference scattering angle 6 = 7r/2, while here it is much 
larger. It is easy to understand this result considering the structure of the MHV amplitudes in the 
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Figure 3: Differential distributions for the minimum gluon transverse momentum. Exact 
result vs. SPHEL. 

two cases. In the gluon-only case, eq. (Effi ), the numerators are dominated by the term (P1P2) 4 ~ <s 4 , 
while in the qq case, eq. (|31|), only terms proportional to t-channel momentum exchange appear. 
With the increase in the number of final-state particles, the average momentum exchanged in the 
t-channels becomes smaller, while s stays the same, and the ratio s/ (t) 4 becomes very large. 

We also considered differential cross-sections. In Fig. 3, we show the distribution of the mini- 
mum gluon transverse momentum for both processes. A good agreement between the exact and 
the SPHEL result is observed. Considering the large uncertainties present in the overall normal- 
ization of multi-parton tree-level processes (e.g. those induced by changes in the choice of the 
renormalization and factorization scales), the approximated evaluation can provide in many cases 
a sufficient description of the final-state distributions. 

In Figs. 4 and 5 we plot the distributions for the maximum gluon transverse momentum and 
the maximum two-gluon invariant mass. Again there is a reasonable agreement between the exact 
and the approximated result, in particular in the case of the gluon-only process. 

Before closing this Section, we present some technical details to illustrate the performance of 
the algorithms used to produce our results. As an example, consider the process g g — ► 8 g. As 
shown in Table 1, the total number of Feynman diagrams contributing to this process is 10,525,900. 
The plots in this work have been obtained from the evaluation of the matrix elements for 1.9 x 10 6 
Monte Carlo points passing the selection cuts given in eq. ( |32"D (out of 1.8 x 10 8 points selected 
by the phase-space generator). The efficiency for the generation of the non zero weights, defined 
as the average weight divided by the maximum weight, was about 1 x 10~ 4 . The computational 
time for producing 100 events that pass the cuts, with a 200 MHz processor 0, and working in 
double precision throughout, was about 91 seconds for the exact matrix element and 2.3 seconds 
using the SPHEL approximation. In this last case, the dominant component of the CPU budget is 

11 All of the following numbers are reduced by a factor of 4 when using a DEC Alpha station. 
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Figure 4: Differential distributions for the maximun gluon transverse momentum. Exact 
result vs. SPHEL. 
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the search for phase-space configurations passing the selection cuts. 



5 Conclusions 

We presented in this work an algorithm to evaluate the exact, tree-level matrix elements for 
multi-parton processes in QCD. This technique, based on the algorithm ALPHA, has been tested 
for processes such as gg — > n gluons and qq —> n gluons, with n up to 9. We discussed how 
the summation over colour configurations allows the construction of parton-level event generators 
suitable to interfacing with a parton-shower evolution including the effects of colour-coherence. 
This will eventually lead to a fully exclusive, hadron-level description of multi-jet final states, 
accurately incorporating the dynamics of large jet-jet separation angles. 

While we confined our presentation to the case of hadroproduction, our program can be used 
without modifications for the evaluation of multi-parton production in e + e~ collisions, since the 
relevant SU{2) x U{\) Lagrangian is already included in the code. Likewise, associated hadropro- 
duction of gauge bosons and QCD partons is a straightforward application of our code. 

We gave some explicit numerical results, considering parton-level cross-sections and distribu- 
tions for the processes gg — > 8 gluons and qq — > 8 gluons. We verified that some standard simple 
approximations to the multi-parton cross-sections provide a good description of the exact result, 
for both rates and distributions. The existence of the exact results allows now to extend the checks 
on these approximations to values of n larger than ever before. We expect that large improvements 
can be obtained in the numerical efficiency of the algorithm, and that cross- sect ions for up to 10 
final-state partons should become feasible. 

Furthermore, when n becomes so large that the evaluation of the exact amplitudes is nu- 
merically too slow for the generation of large samples of events, one can nevertheless use the 
approximated calculations to generate the samples, and lower statistics runs to assess the good- 
ness of the approximation. In this respect, the numerical efficiency of the phase-space generation 
in presence of kinematic cuts will need to be improved as well. 
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